Ursus Americanus - Forest Cover Analysis

# load and view data
load("bc_bear_occurrences.Rda")
# str(occ_data)
load("BC_Covariates.Rda")
summary(DATA)

##            Length Class           Mode
## Window      1     SpatialPolygons S4  
## Elevation  10     im              list
## Forest     10     im              list
## HFI        10     im              list
## Dist_Water 10     im              list

Visualize bear occurrences in BC

# extract location columns
bears_loc <- occ_data[, c("decimalLongitude", "decimalLatitude", "month", "year")]
bears_loc_filtered <- subset(bears_loc, year %in% c(2020, 2021, 2022, 2023, 2024))

# create sf object
bears_sf <- st_as_sf(
  bears_loc_filtered,
  coords = c("decimalLongitude", "decimalLatitude"),
  crs = 4326  # WGS84 (longitude/latitude)
)

# BC Albers projection
bears_sf_proj <- st_transform(bears_sf, crs = 3005)

# extract x, y coordinates
coords <- st_coordinates(bears_sf_proj)

# extract BC window
window_sf <- st_as_sf(DATA$Window) # convert SpatialPolygons to Simple Features (sf)
window_proj <- st_transform(window_sf, crs = 3005) # ensure same CRS
window <- as.owin(window_proj) # convert to owin using sf object

# create ppp object
bears_ppp <- ppp(
  x = coords[,1],
  y = coords[,2],
  window = window
)

## Warning: 187 points were rejected as lying outside the specified window

## Warning: data contain duplicated points

plot(bears_ppp, pch = 21, main = "Black bear occurrences in BC, 2020-2024")

## Warning in plot.ppp(bears_ppp, pch = 21, main = "Black bear occurrences in BC,
## 2020-2024"): 187 illegal points also plotted

Compare seasonal distributions

# Define the seasons
seasons <- list(
  winter = c(12, 1, 2),
  spring = c(3, 4, 5),
  summer = c(6, 7, 8),
  autumn = c(9, 10, 11)
)

# Create an empty list to store the ppp objects for each season
ppp_list <- list()

# Create an empty list to store the filtered data for each season
bears_sf_list <- list()

# Loop over seasons
for (i in 1:length(seasons)) {
  
  # Filter for each season
  season_name <- names(seasons)[i]
  season_months <- seasons[[i]]
  
  # Filter bears_loc_filtered by the season's months
  bears_loc_season <- bears_loc_filtered[bears_loc_filtered$month %in% season_months, ]
  
  # Print the number of observations for this season
  cat(season_name, ": ", nrow(bears_loc_season), " observations\n", sep="")
  
  # Create sf object for the filtered season
  bears_sf_season <- st_as_sf(
    bears_loc_season,
    coords = c("decimalLongitude", "decimalLatitude"),
    crs = 4326  # WGS84 (longitude/latitude)
  )
  
  # Transform to BC Albers projection
  bears_sf_season_proj <- st_transform(bears_sf_season, crs = 3005)
  
  # Store sf object in the list
  bears_sf_list[[season_name]] <- bears_sf_season_proj
  
  # Extract x, y coordinates
  coords <- st_coordinates(bears_sf_season_proj)
  
  # Create ppp object for each season
  suppressWarnings(
    bears_ppp_season <- ppp(
      x = coords[, 1],
      y = coords[, 2],
      window = window
    )
  )
  
  # Store ppp object in the list
  ppp_list[[season_name]] <- bears_ppp_season
}

## winter: 76 observations
## spring: 808 observations
## summer: 1862 observations
## autumn: 841 observations

# Plot the four maps in a 2x2 grid
par(mfrow=c(2,2), mar=c(1,1,1,1))  # Set up a 2x2 plotting window

# Loop over ppp objects and plot them
for (i in 1:length(ppp_list)) {
  season_name <- names(ppp_list)[i]
  suppressWarnings(
    plot(ppp_list[[season_name]], main = season_name, pch = 21)
  )
}

Looking at covariates

par(mfrow=c(2,2), mar=c(1,1,1,1))  # Set up a 2x2 plotting window
plot(DATA$Elevation)
plot(DATA$Forest)
plot(DATA$HFI)
plot(DATA$Dist_Water)

Forest cover analysis

forest <- DATA$Forest

# intensity as a function of forest cover
par(mfrow=c(2,2), mar=c(2,2,2,2))
for (i in 1:length(ppp_list)) {
  season_name <- names(ppp_list)[i]
  rho <- rhohat(ppp_list[[season_name]], forest)
  plot(rho, xlim=c(0, max(forest)), main = paste(season_name))
}

# forest cover for bears locations

im2raster <- function(im_obj) {
  m <- as.matrix(im_obj) # image to matrix
  r <- raster(m) # create raster object
  # range and extension
  ext <- c(im_obj$xrange[1] - im_obj$xstep/2,
           im_obj$xrange[2] + im_obj$xstep/2,
           im_obj$yrange[1] - im_obj$ystep/2,
           im_obj$yrange[2] + im_obj$ystep/2)
  extent(r) <- ext
  return(r)
}

# extract forest cover values at occurrence locations
forest_raster <- im2raster(forest)
bears_sf_proj$forest_value <- extract(forest_raster, bears_sf_proj)

KDE to estimate the concentration of bear occurrences based on forest cover values.

# KDE for Bear Locations
kde_b <- density(bears_sf_proj$forest_value, na.rm = TRUE, bw = "SJ-dpi")
summary(bears_sf_proj$forest_value)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    0.00   25.44   52.36   53.18   81.17  100.00    1503

plot(kde_b,
     main = "Kernel Density Estimate of Forest Cover at Bear Locations",
     xlab = "Elevation (meters)",
     ylab = "Density",
     col = "#2C6B2F", 
     lwd = 2)  

These results suggest that bears are found in a mix of forested and non-forested areas, but on average, they’re in areas with moderate to high forest cover.

Modelling

Fitting bear location as a quadratic function of forest cover

# NA with mean forest cover
mean_forest <- mean(as.vector(as.matrix(forest)), na.rm = TRUE)
forest_clean <- eval.im(ifelse(is.na(forest), mean_forest, forest))

# polynomial model
forest_model <- ppm(bears_ppp ~ forest + I(forest^2), covariates = list(forest = forest_clean))
print(forest_model)

## Nonstationary Poisson process
## Fitted to point pattern dataset 'bears_ppp'
## 
## Log intensity:  ~forest + I(forest^2)
## 
## Fitted trend coefficients:
##   (Intercept)        forest   I(forest^2) 
## -2.034720e+01  6.209799e-02 -6.168772e-04 
## 
##                  Estimate         S.E.       CI95.lo       CI95.hi Ztest
## (Intercept) -2.034720e+01 4.834276e-02 -2.044195e+01 -2.025245e+01   ***
## forest       6.209799e-02 2.093576e-03  5.799466e-02  6.620132e-02   ***
## I(forest^2) -6.168772e-04 2.010334e-05 -6.562791e-04 -5.774754e-04   ***
##                   Zval
## (Intercept) -420.89447
## forest        29.66120
## I(forest^2)  -30.68531

# predicts and plots intensity
pred_intensity_forest <- predict(forest_model)
plot(pred_intensity_forest, 
     main = "Predicted Intensity Based on Forest Cover")
points(bears_ppp, 
       pch = 20, 
       col = "red")

The positive coefficient on the linear term indicates that bear occurrence intensity initially increases with forest cover. However, the negative coefficient on the quadratic term means that this relationship has a peak, after which further increases in forest cover are associated with a decrease in expected intensity. This suggests that bears are most likely to be found in areas with moderate forest cover, and less likely to occur in areas with either very low or very high forest density.

This pattern may reflect ecological preferences, where bears favor mixed habitats over open or densely forested areas. It’s also possible that human detection rates are lower in areas of dense vegetation, where visibility is reduced.

Computing and plotting residuals

forest_residuals <- residuals(forest_model, type = "pearson")
forest_residuals$v[!is.finite(forest_residuals$v)] <- NA # removing any non-finite residuals
plot(forest_residuals, 
     main = "Residuals - Forest Model", 
     na.col = "transparent")

Seasonal analysis

Fitting a model for each season

models_seasonal <- list()

# model for each season (individually)
for (season in names(ppp_list)) {
  models_seasonal[[season]] <- ppm(ppp_list[[season]], ~ forest + I(forest^2), covariates = list(forest = forest_clean))
  # cat("Model for", season, ":\n")
  # print(models_seasonal[[season]])
  plot(predict(models_seasonal[[season]]), main = paste("Intensity ~ Forest Cover -", season))
}

## Warning: glm.fit: algorithm did not converge

Winter: Winter: Intercept: -23.78; Forest Coefficient: 0.0582; Quadratic Coefficient: -0.000672.

Bear occurrences increase with forest cover up to a point, then decline. The smaller linear coefficient suggests a slower rise in intensity compared to other seasons. The model did not converge, so caution is needed in interpreting the results.

Spring: Intercept: -21.96; Forest Coefficient: 0.0657; Quadratic Coefficient: -0.000632.

The relationship is similar to winter but slightly stronger, meaning a more noticeable increase and peak in moderate forest areas.

Summer: Intercept: -21.05; Forest Coefficient: 0.0621; Quadratic Coefficient: -0.000602.

Summer shows a comparable trend with a slightly weaker curvature, suggesting bears may tolerate higher forest density before intensity starts to drop.

Autumn: Intercept: -21.70; Forest Coefficient: 0.0678; Quadratic Coefficient: -0.000732.

Autumn has the sharpest curve (greatest linear and quadratic coefficient magnitudes), indicating a more distinct peak and a stronger decline in high-density forest areas.

Overall, all seasons display a consistent non-linear relationship where bear occurrence intensity increases with forest cover up to a moderate level, then declines. While the shape of the curve remains similar, the strength and location of the peak vary slightly by season, indicating subtle shifts in habitat preference throughout the year.

Fitting bear location as a quadratic function of forest cover with season as a factor

# Combine the data from each season into a single data frame
bears_all <- do.call(rbind, lapply(names(bears_sf_list), function(season) {
  sf_season <- bears_sf_list[[season]]
  coords <- st_coordinates(sf_season)
  forest_vals <- extract(forest_raster, sf_season)
  data.frame(x = coords[, 1], y = coords[, 2],
             forest_value = forest_vals, season = season)
}))

# Convert "season" to a factor
bears_all$season <- as.factor(bears_all$season)

# Create a ppp object for the entire dataset using the defined study window.
bears_ppp_all <- ppp(
  x = bears_all$x,
  y = bears_all$y,
  window = window,
  marks = bears_all$season
)

## Warning: 187 points were rejected as lying outside the specified window

## Warning: data contain duplicated points

# Check how many points remain in the ppp object
cat("Number of points in bears_ppp_all:", npoints(bears_ppp_all), "\n")

## Number of points in bears_ppp_all: 3400

# Retrieve the marks from the ppp object (these correspond to the points kept after rejection)
current_marks <- marks(bears_ppp_all)

# Convert the marks into a data frame with the correct number of rows
marks(bears_ppp_all) <- data.frame(season = current_marks)

# Fit the combined model including the interaction between HFI and season.
mod_combined <- ppm(bears_ppp_all ~ forest + I(forest^2) + marks, 
                    covariates = list(forest = forest_clean))

# Display the model summary
mod_combined

## Nonstationary multitype Poisson process
## Fitted to point pattern dataset 'bears_ppp_all'
## 
## Possible marks: 'autumn', 'spring', 'summer' and 'winter'
## 
## Log intensity:  ~forest + I(forest^2) + marks
## 
## Fitted trend coefficients:
##   (Intercept)        forest   I(forest^2)   marksspring   markssummer 
## -2.179537e+01  6.209799e-02 -6.168772e-04 -3.957121e-02  7.891398e-01 
##   markswinter 
## -2.379296e+00 
## 
##                  Estimate         S.E.       CI95.lo       CI95.hi Ztest
## (Intercept) -2.179537e+01 5.739746e-02 -2.190787e+01 -2.168287e+01   ***
## forest       6.209799e-02 2.093576e-03  5.799466e-02  6.620132e-02   ***
## I(forest^2) -6.168772e-04 2.010334e-05 -6.562791e-04 -5.774754e-04   ***
## marksspring -3.957121e-02 5.053363e-02 -1.386153e-01  5.947288e-02      
## markssummer  7.891398e-01 4.266227e-02  7.055233e-01  8.727563e-01   ***
## markswinter -2.379296e+00 1.215116e-01 -2.617454e+00 -2.141137e+00   ***
##                     Zval
## (Intercept) -379.7270268
## forest        29.6612018
## I(forest^2)  -30.6853099
## marksspring   -0.7830669
## markssummer   18.4973701
## markswinter  -19.5808065

In general, as forest cover increases, the log intensity of bear occurrence initially increases. The negative quadratic term indicates that the intensity eventually decreases after a peak, suggesting a preference for moderate forest cover.

The magnitude of the spring coefficient is small compared to autumn (baseline), suggesting that after adjusting for forest cover, spring does not differ significantly from the autumn in terms of intensity of bear occurrences.

Summer has the highest relative intensity among seasons, meaning bear occurrences are much more likely in summer compared to autumn.

Winter has a much lower intensity than autumn, indicating far fewer bear sightings or activity.

Performing Quadrat test for each season

# quadrat tests
for (s in names(ppp_list)) {
  qt <- quadrat.test(ppp_list[[s]], nx = 5, ny = 5)
  print(qt)
  plot(qt, main = paste("Quadrat test -", s))
}

## Warning: Some expected counts are small; chi^2 approximation may be inaccurate

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  ppp_list[[s]]
## X2 = 561.39, df = 20, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 21 tiles (irregular windows)

## Warning: Some expected counts are small; chi^2 approximation may be inaccurate

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  ppp_list[[s]]
## X2 = 2202.7, df = 20, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 21 tiles (irregular windows)

## Warning: Some expected counts are small; chi^2 approximation may be inaccurate

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  ppp_list[[s]]
## X2 = 3767, df = 20, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 21 tiles (irregular windows)

## Warning: Some expected counts are small; chi^2 approximation may be inaccurate

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  ppp_list[[s]]
## X2 = 2720.5, df = 20, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 21 tiles (irregular windows)

All p-values are small (< 2.2e-16) indicating significant deviation from CSR in each season.